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Contact and Quasi-Static Impact of a Dissipationless Mechanical Model 
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Collisions and contacts of elastic materials are numerically and theoretically investigated. 
Using a two-dimensional spring-mass model with defect particles under the free boundary 
condition, we reproduce the Hertzian contact theory at equilibrium and the quasi-static theory 
for low speed impacts. 
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1. INTRODUCTION 

The origin of irreversibility from a reversible mechan- 
ical model is one of the most fundamental subjects in 
non-equilibrium statistical mechanics. A typical exam- 
ple can be seen in transport phenomena in nonlinear 
spring models. 1,2 Although most of models discuss the 
transport processes under the influence of thermostats, 
we still do not understand the mechanism to reach an 
equilibrium state based on a purely mechanical model. 

The aim of this paper is to derive macroscopic laws 
of elastic materials based on a microscopic lattice model. 
Here we discuss contacts and impacts between elastic ma- 
terials. The former is related to the origin of the second 
law of thermodynamics from a purely mechanical model. 

For contacts of elastic materials, we believe that the 
contacts between elastic bodies can be described by 
the Hertzian contact theory. 3-5 The two-dimensional 
Hertzian contact theory gives us the relation between 
the deformation of a disk 8 and the compressive force P 
as 
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(1) 



where R is the radius of the undeformed disk. 6 Here, 
E* is the reduced Young's modulus, E* = E/(l - v 2 ), 
where E and v are Young's modulus and Poisson's ratio, 
respectively. 

For the low speed head-on collisions of elastic materi- 
als, the relation between the restitution coefficient e and 
the impact speed v is described by the quasi-static the- 
ory of low speed impacts. 7,8 The quasi-static theory is 
an extension of the Hertzian contact theory to include 
the internal viscosity of materials. By solving the equa- 
tion of motion for the deformation with adequate initial 
conditions and calculating the rebound speed, we can 
obtain the relation between the impact speed v and the 
restitution coefficient e. 

The restitution coefficient e depends also on the in- 
cident angle. 9,10 We have recently carried out a two- 
dimensional simulation of oblique impacts using a dis- 
sipationless elastic model based on the mass-spring 
model 11, 12 to reproduce and explain the previous experi- 
mental result by Louge and Adams. 9 Although the model 
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can reproduce similar results to the experimental results, 
one of the apparent defects in the model is that the sys- 
tem cannot reach an equilibrium state described by the 
Hertzian contact theory. 3-5 When we put the disk on the 
wall under the influence of the gravity, the release of ini- 
tial potential energy is recurrent and the system keeps os- 
cillations. We also have expected that e can be described 
by the quasi-static theory for low speed head-on colli- 
sion of the disk with a potential wall, but we have found 
a discrepancy between the numerical result and the pre- 
diction by the quasi-static theory. 6, 13 Refining the model 
to overcome those defects in a dissipationless model will 
be helpful to understanding the origin of irreversibility 
and the second law of thermodynamics from a mechani- 
cal point of view. 

In this paper, we propose a microscopic dissipation- 
less model to reproduce the Hertzian contact at equilib- 
rium and the quasi-static theory for low-speed impacts. 
These results enable us to understand the origin of irre- 
versibility for finite degrees of freedom from a reversible 
mechanical model. The construction of this paper is as 
follows. In the next section we introduce our model. In 
Sec. 3 and 4 we show the results of our simulation of con- 
tact problems and low-speed impacts, respectively, and 
briefly summarize our results. Appendix A is devoted to 
the calculation of elastic moduli of our model. 

2. Model 

Let us introduce our model(Fig.l). It is basically the 
same as our previous model which obeys Hamilton equa- 
tion. 11,1 The disk consists of 1099 mass points while 
the wall consists of 1269 mass points. The two corner 
points of the bottom of the wall are fixed. The surfaces 
of both the disk and the wall are initially flat. The in- 
teraction between the disk and the wall is introduced as 
follows. Each mass point i on the lower half boundary of 
the disk receives the force, F(lg ) = aVo exp(— a?i I ' l )ni 4 ' ) , 

(i) 

where l s is the distance between i-th surface mass point 
of the disk and the nearest surface spring of the wall, 
a = 500/ R, Vq — amc 2 R/2, m is the mass of each mass 
point i, cis the one-dimensional speed of sound, and n s 
is the unit vector normal to the connection between two 
surface mass points of the wall. 11, 14 Thus, the dynamical 
equation of motion for each mass point i of the lower half 



2 J. Phys. Soc. Jpn. 



Full Paper 



Author Name 




4R 



Fig. 1. The clastic disk and wall consisted of random lattice. We 
set y = at the surface of the wall. In this figure, the height 
and the width of the wall are respectively R and 4R with the 
undeformed radius of the disk R. 



boundary of the disk is described by 
d 2 r 



ni- 



di 2 



fc 6 x?.}+e(/ tfc -ZW)aV exp(-a/« 

(2) 

where is the position of z-th mass point, t is the time, 
Ni is the number of mass points connected to i-th mass 
point, Xij is the relative deformation vector of the spring 
from the natural length between i-th and j-th connected 
mass points, k a and fc& = k a x 10 _3 /i? 2 are the spring 
constants. Here 0(x) is the step function, i.e. Q(x) = 1 
for x > and Q(x) — for x < 0, and the threshold 
length l t h is the average of the natural lengths of the 
springs of the disk. For internal mass points, the last 
term of the right hand side of eq.(2) is omitted. In most 
of our simulations, we adopt k a = = 1.0 x mc 2 / R 2 
for the disk and k a — k[ w ^ = 1.0 x I0 2 mc 2 /R 2 for the 
wall. Numerical integration of eq. (2) is performed by the 
fourth order symplectic integrator with the time step 
dt = lO~ 3 R/c. 

There are two main differences between our new model 
and the previous model. Here we introduce some defect 
particles and adopt the free boundary condition in which 
the kinetic energy can escape from the boundary. The de- 
fect particles are introduced by choosing mass point i at 
random and eliminating Ni — 1 connecting bonds among 
Ni connecting bonds to i-th mass point. We have intro- 
duced 10 defects for each body. The reason to choose 10 
defects will be discussed later. With the introduction of 
these defects, we expect that the motion of defects be- 
comes irregular and the vibrating wave can be localized 
without spreading. This irregular motion may create the 
irreversibility of time evolution of mass points. In ad- 
dition, we adopt the free boundary condition for both 
sides and the bottom of the wall in contrast to the re- 
flective boundary condition in our previous model. We 
regard the wall as a part of a large system. When we 
put a mechanical perturbation in the wall such as a con- 
tact or an impact, the effect propagates as elastic waves 
and goes out from the boundary because our system is 
a part of a larger system. In this situation, J • > 



should be satisfied where n?, and J are respectively the 
unit normal vector on the boundary and the energy flux 
J = X)i{(l/2) TOV i + £i}vi with the potential energy e*. 
In order to realize the condition, at each time step of 
numerical integration, we set = if Vj is directed to 
inner region of the wall. Because the reflected wave does 
not come into the system, the total energy of our system 
is not conserved. 

Here we briefly comment on the boundary condition 
we adopt. Our model does not conserve the total en- 
ergy because there is the outgoing energy flux from the 
system. We can put thermostats, e.g. Nose-Hoover ther- 
mostats, 15 on the boundary to keep the energy conserva- 
tion law. However, once we introduce the thermostat in 
the system, the phase volume is not conserved, and there 
is the entropy production. 16 In addition, we believe that 
the thermalization from the thermostat does not play an 
important role in macroscopic elastic materials. 13 Thus, 
we adopt the model without the conservation of the en- 
ergy. In the next section we will show that the introduc- 
tion of the thermostat does not make much difference in 
our results. 

3. Simulation of Contact Problems 

In this section, we show the results of our simulation 
of elastic contacts. At first, we show the results of the 
simulation for the contact problem. In this simulation, 
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Fig. 2. The relation between the deformation and the external 
force. Cross points with error bars are the numerical results with 
the initial temperature T = 0. Plus points are the numerical 
results with the initial temperature T = 10 — 7 Mc 2 in the simu- 
lation with Nose-Hoover thermostats. The solid line is the pre- 
diction by the Hertzian contact theory with v = 0.336. 



we introduce the wall whose height and width arc AR 
and R, respectively. We put the disk in the external field 
P ranging from P = 5.77 x 10- 3 ttRE* to P = 1.27 x 
1Q- 2 iiRE*. Thus, in this case, we add the term —(P/N)y 
in the right hand side of eq.(2), where N — 1099 and y 
is the unit vector in y direction. The initial oscillation of 
the center of mass of the disk relaxes to reach a stable 
oscillation. After the relaxation of the center of mass, 
we calculate the deformation of the disk, 8 = R — Rd, 
where Rd is the deformed radius of the disk which is the 
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distance from the contact patch to the center of mass of 
the disk. 

Figure 2 is the relation between 5/R and P/ttRE*. 
Cross points are the averaged results of 10 disks with dif- 
ferent configuration of mass points, and the error bars are 
the standard deviations. The solid line in Fig. 2 is eq.(l) 
with v — 0.336 calculated from the two-dimensional the- 
ory of elasticity based on the assumption of an isotropic 
disk. Details of the calculation of elastic constants of the 
model are shown in Appendix A. Our simulation data 
are well reproduced by the theoretical curve which does 
not include any fitting parameters. Thus, we conclude 
that our model can produce an equilibrium state at the 
contact. 

Let us discuss the effect of thermalization from ther- 
mostats on our results. We have carried out another sim- 
ulation, in which the wall has fixed boundaries at the 
bottom and the both side ends, and the mass points 
connected to the fixed boundaries obey the equation of 
motion of particles connected to the Nose-Hoover ther- 
mostat. 15 We have arranged that the temperature of the 
wall becomes T = 10~ 7 Mc 2 . Here we define the temper- 
ature of the thermostat by the variance of the Maxwell- 
Boltzmann distribution which we give as the initial veloc- 
ity distribution of the thermostat particles. Plus points 
in Fig. 2 are the averaged results of 10 disks with differ- 
ent configuration of mass points in the simulation with 
the Nose-Hoover thermostat. This result shows that the 
introduction of thermostats does not affect the relation 
between the compressive force and the deformation so 
much. From this result, we conclude that the thermliza- 
tion from the thermostat does not play an important role 
in our simulation of contact problems. 

Introducing defect particles in the model plays an im- 
portant role in the relaxation of internal vibration in 
the contact problem. To characterize the relaxation pro- 
cess to an equilibrium state, we investigated the time 
evolution of the velocity distribution function (VDF) 
f{v x ,v y ,t), where v x and v y are velocity components 
of mass points of the disk, and the Shannon entropy 
when the compressive force P = 8.75 x 10~ 3 7r RE*. To 
obtain the velocity distribution function f(v x ,v y ,t) = 
J2xJ2 y f(x,y,v x ,v y ,t), we calculate f(x,y,v x ,v v ,t) by 
averaging the data from t — 6R/c to t. 

Figure 3 shows the time evolution of VDF f(v x , v yi t). 
The distribution function near the initial stage (Fig. 
3(a)), f(v x ,v y7 t = 6R/c), has three peaks along the 
axis v x = 0. The peaks may be attributed to the coex- 
istence of the reflective motion and the compressive mo- 
tion of mass points in the disk. Figure 3(b) shows that 
f{v X} v y ,t max = 60R/c) has a Gaussian form. In fact, 
f{v Xl 0, t max ) can be well fitted by the Gaussian whose 
variance is 0.02. 

In addition, we have investigated the time evolu- 
tion of the Shannon entropy. Figure 4 shows the 
time evolution of the entropy defined by S(t) = 
- T, x T, y Et, x E„ f( x , V, v x ,v y , t) In f(x, y, v x ,v y , t) for 
several numbers of defects. We have investigated the de- 
pendence of the number of defect particles on the en- 
tropy. In Fig. 4, we change the number of defects in the 
disk from to 15, and investigate the time evolution of 
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Fig. 3. The velocity distribution of mass points of the disk at (a) 
t = to and (b) t = t max ■ 

the entropy for each disk. In the case of the disks whose 
number of defects are from 8 to 10, after t — 30R/c, the 
time evolution of the entropy is described as an univer- 
sal curve which shows monotonic increase and reaches 
the maximum value 13.5. In the case of other disks, the 
entropy does not show the monotonic increase. The mix- 
ing effect is not enough for the smaller number of defects 
while a global oscillation is excited because of the weak- 
ness of the disk structure for too large number of defects. 

We have also calculated the averaged kinetic energy 
per one defect particle (Fig. 5). The time averaged ki- 
netic energy of one defect is not so large as that of the 
other particles. However, the kinetic energy for one defect 
changes at random with larger amplitude. This tendency 
originates from irregular motion of the defect particles, 
which enables the oscillation modes to mix and to reach 
an equilibrium state. In realistic systems, such a phe- 
nomenon can be seen in phonon scattering by defects 
and impurities of solids, the scattering which realizes a 
thermal equilibrium state. We conclude that the irregu- 
lar motion of the defect particles generates irreversibility 
and makes the model reach an equilibrium state, which 
reproduces the Hertzian contact of our model. 

4. Simulation of Low-speed Impact 

In this section, we show the results of the simulation of 
low speed impact. Here we arrange that the height and 
the width of the wall are respectively R and AR. From 
the simulation of the head-on collision of the disk with 
the wall, we calculate e for each initial speed. The initial 
speed is ranged from w = 1.0x 10~ 3 c to v = 1.0 x 10~ 2 c 
without the external field. 

Figure 6 (a) is the relation between the restitution 
coefficient e and the impact speed v/c ranging from 
v/c = 0.001 to v/c — 0.02 while Fig. 6 (b) is the results in 
the low colliding speed. The cross points are the averaged 
results of 10 disks with different configuration of mass 
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Fig. 4. The time evolutions of the Shannon entropy S(t) when 
the number of defects are (a) from 8 to 10 and (b) and 15. X- 
axis is the scaled time while y-axis is the Shannon entropy S(t) 
calculated from the distribution function of internal mass points 
of the disk. 
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Fig. 5. The time evolution of the kinetic energy per one particle 
of bound particles and defect particles of the disk. 



points, and the error bars are the standard deviations. 
As can be seen, e decreases slightly with increasing collid- 
ing velocity and e cannot reach 1 as impact velocity de- 
creases. We compare this result with the two-dimensional 
quasi-static theory of low speed impact. 7 ' 8,17 In the two- 
dimensional quasi-static theory, 18 the dynamical equa- 
tion of motion for the deformation of the disk may be 
written as 
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Fig. 6. The relation between the colliding speed v and the resti- 
tution coefficient e. The solid line indicates the prediction by the 
quasi-static theory. The cross points with error bars are numer- 
ical results in the range (a) from v/c = 0.001 to v/c = 0.02 and 
(b) from v/c = 0.001 to v/c = 0.007. 



where M is the mass of the disk and To is the time scale 
of dissipation. P is obtained as a function of 8 by numeri- 
cally solving eq.(l) for P. The last term in the right hand 
side is the dissipative force which is proportional to the 
velocity of the deformation. 7 ' 8 By introducing To as a fit- 
ting parameter, we solve eqs.(l) and (3) numerically with 
the initial conditions (5 = and dS/dt = v to obtain the 
relation between e and v. The solid line in Fig. 6 (a) and 
(b) is the numerical result of eq.(3) with To = 0.011i?/c. 
Figure 6 (b) shows that our simulation data are well re- 
produced by the quasi-static theory in the low speed re- 
gion ranging from v/c = 0.001 to v/c = 0.007. However, 
Fig. 6 (a) shows that the quasi-static theory is no longer 
valid in the high speed region larger than v/c= 0.007 in 
our simulation. The excitation of various internal modes 
originated from high speed impact causes such the dis- 
crepancy. 

5. Concluding Remarks 

In summary, we reproduce Hertzian contact mecha- 
nism without introduction of any explicit dissipations. 
Our numerical results for low-speed impact are also 
consistent with the two-dimensional quasi-static theory. 
Through our investigation, we expect that we can ob- 
tain the further understanding on the origin of the irre- 
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versibility from the reversible kinetic model. 

Our future tasks are (i) deriving the characteristic time 
r in quasi-static impact from the microscopic mechan- 
ical model, (ii) extending our analysis to three dimen- 
sional problem, (iii) constructing a theory of impact for 
non quasi-static region, (iv) investigating the effect of the 
density distribution of colliding materials on their impact 
behavior, and (v) to check the effect of gravity for quasi- 
static impacts. For the last point, our preliminary results 
suggest that the restitution coefficient e depends on the 
value of the gravitational acceleration. 
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Appendix: Calculation of Elastic Moduli 

In this appendix, we show how to calculate elastic 
moduli, such as Young's modulus E and Poisson's ratio 
v, of our model. In a two-dimensional isotropic medium, 
the energy density E p becomes 

E p = 2e 2 (\ + f i) (A-l) 

in an isotropic compression while 

E p = l/2^ie 2 (A-2) 

in a simple shear, where e is the strain, and A and v are 
Lame's constants. 

We calculate the energy densities of the disk in our 
model by adding an isotropic compression and a simple 
shear with the strain e. For an isotropic compression, by 
adding the small strain e to each spring of the disk with 
the natural length dj , we calculate the energy density of 
the disk as 

^( dj ef + ^(djeA (A-3) 



3 

j 



where A is the of the disk. We neglect the fourth- 
order term. From eqs.(A-l) and (A-3), on the assumption 
that our model is isotropic, we can calculate A + [i as 
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(A-4) 



For a simple shear, by displacing the position of each 
mass point (a;,, yi) as (x t + ey i} yt) with the small strain 
e, we calculate the energy density of the disk as 



E„ 
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(A-5) 



length of the j-th spring. Sdj is calculated as 

Sdj = yf {xi - x\ + e(yl - j^')} 2 + {yi - y^') 2 - dj 



{<- x{){y 3 a -yl), 



(A-6) 



where (x J a ,y 3 a ) and (x J b ,yl) are the initial positions of 
both ends of the j-th spring. Thus, cq.(A-5) becomes 

E P ce 2 ^j2^<-<) 2 (yi-yl) 2 - ( A - 7 ) 

From eqs.(A-2) and (A-7), we can calculate \x as 
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2A ^ d 2 

3 3 
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(A-8) 



Technically, we input the data for initial configura- 
tion of the connecting bonds of the disk into MATHE- 
MATICA and put perturbations, such as the isotropic 
compression and the simple shear, to calculate the en- 
ergy densities. Calculating A and /i from cqs.(A-4) and 
(A-8) to substitute them into the two-dimensional rela- 
tion, E = 4/i(A + A*)/(2/i + A) and v = A/(2/i + A), we 
calculate Young's modulus E and Poisson's ratio v as 
E = 0.773fc a and v = 0.336, respectively. 
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where Sdj is the relative displacement from the natural 



